function [likelihood_contribution]=likelihood_contrib(parameter,ability_metric,LowCost_Ability,HighCost_Ability,V_l,V_h,outcome,tpORtp3h,Cpm25,parameter_fixed,LowCost_Error,HighCost_Error)
% Compared to likelihood_obj, outputs Nx1 vector of likelihood contributions

%N=size(LowCost_Ability,1);

% Parameters that are fixed
%cutoff_tp=parameter_fixed(1);
%cutoff_pm=parameter_fixed(2);
%delta_0  =parameter_fixed(3);

% Parameters that are estimated 
k       =parameter(1);
%delta_tp=parameter(3);
%delta_pm=parameter(4);
                
% Delta changes separately in the log of (1 +) the temperature or PM2.5 difference relative to a cutoff, i.e., temperature of 27 -> log (1 + 27 - cutoff) = log (1 + 27 - 27) = log 1 = 0; 28->.7 of delta; 43->2.8   
%delta = repmat(delta_0,N,1) + delta_tp * log( max( (tpORtp3h - repmat(cutoff_tp,N,1)) , 0) + ones(N,1)) + delta_pm * log( max( (Cpm25 - repmat(cutoff_pm,N,1)) , 0) + ones(N,1));


% Marginal cost (mc) function:
% inputs: ability metric, cost parameter, LowCost_Ability, HighCost_Ability
% outputs: c_l, c_h
[c_l,c_h]=mc(ability_metric,parameter,LowCost_Ability,HighCost_Ability,LowCost_Error,HighCost_Error);

% transition_prob_oneshot:
% inputs: c_l, c_h, V_l, V_h, k
% outputs: p_l_winoneshot=p_l_win3, p_h_winoneshot=p_h_win3
[p_lwin3,p_llose3]=transition_prob_oneshot(c_l,c_h,V_l,V_h,k);

% transition_prob_stage2:
% inputs: c_l, c_h, V_l, V_h, temperature, PM2.5, k, parameters of the delta function, calls on transition_prob_oneshot
% outputs: Continuation values (in stage 3) conditional on stage-2 outcome (and stage-1 history), stage-2 transition probabilities  
[~,~,~,~,~,~,~,~,p_lwin2_lwon1,p_llose2_lwon1,p_lwin2_llost1,p_llose2_llost1]=transition_prob_stage2(parameter,tpORtp3h,Cpm25,parameter_fixed,c_l,c_h,V_l,V_h);

% transition_prob_stage1:
% inputs: c_l, c_h, V_l, V_h, temperature, PM2.5, k, parameters of the delta function, calls on transition_prob_stage2
% outputs: Continuation values (in stage 2) conditional on stage-1 outcome, stage-1 transition probabilities   
[~,~,~,~,p_lwin1,p_llose1]=transition_prob_stage1(parameter,tpORtp3h,Cpm25,parameter_fixed,c_l,c_h,V_l,V_h);



% Compute probabilities for each match outcome (this can be sped up if necessary subsequently; e.g., for simulated likelihood)
prob_ll =p_lwin1 .*  p_lwin2_lwon1;
prob_lhl=p_lwin1 .* p_llose2_lwon1.* p_lwin3;
prob_lhh=p_lwin1 .* p_llose2_lwon1.*p_llose3;
prob_hll=p_llose1.* p_lwin2_llost1.* p_lwin3;
prob_hlh=p_llose1.* p_lwin2_llost1.*p_llose3;
prob_hh =p_llose1.*p_llose2_llost1;
prob=[prob_ll prob_lhl prob_lhh prob_hll prob_hlh prob_hh];
likelihood_contribution=sum(prob.*outcome,2);
    
    
end